Bistability of an In Vitro Synthetic Autoregulatory Switch 



Pakpoom Subsoontorn\ Jongmin Kim^'^, Erik Winfree^''^''^* 
Departments of 1 Biology, 2 Computation and Neural Systems, 3 Computer Science 

and 4 Bioengineering 
California Institute of Technology 
{pakpoom, jongmin, winfreejOdna. caltech.edu 

January 5, 2011 



Abstract: The construction of synthetic biochemical circuits is an essential step for developing quantitative understanding 
of information processing in natural organisms. Here, we report construction and analysis of an in vitro circuit with 
positive autoregulation that consists of just four synthetic DNA strands and three enzymes, bacteriophage T7 RNA 
polymerase, Escherichia coli ribonuclease (RNase) H, and RNase R. The modularity of the DNA switch template allowed 
a rational design of a synthetic DNA switch regulated by its RNA output acting as a transcription activator. We verified 
that the thermodynamic and kinetic constraints dictated by the sequence design criteria were enough to experimentally 
achieve the intended dynamics: a transcription activator configured to regulate its own production. Although only 
RNase H is necessary to achieve bistability of switch states, RNase R is necessary to maintain stable RNA signal levels and 
to control incomplete degradation products. A simple mathematical model was used to fit ensemble parameters for the 
training set of experimental results and was then directly applied to predict time-courses of switch dynamics and sensitivity 
to parameter variations with reasonable agreement. The positive autoregulation switches can be used to provide constant 
input signals and store outputs of biochemical networks and are potentially useful for chemical control applications. 
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Introduction 

Within a living cell lies an information processing 
system: the genetic circuits, the community of genes 
that regulate one another, and thus allow the cell to 
express the right genes at the right times. Synthetic 
biology provides a new approach to understand design 
principles underlying intricate and dynamic behaviors 
of natural genetic circuits, by building and analyzing 
synthetic circuits which exhibit analogous behaviors; 
this also lays the foundation for future engineering of 
complex chemical and biological systems. With such 
circuits, it is possible to test hypotheses by construc- 
tion, and often synthetic simplicity facilitates quanti- 
tative analysis as well as systematic engineering de- 
sign [imisiiii]. 

For designing and constructing synthetic biochem- 
ical networks, several decision steps are necessary to 
choose the regulatory molecules and the biochemical in- 
frastructure that supports network operation. Protein- 
based synthetic circuits can take advantage of the huge 
diversity of protein structures and functions that al- 
lows a wide range of possible regulatory features O 
[H [71 O [32] ■ Still, from an engineering perspective, it 
remains a challenge to rationally design a new regula- 
tory protein with desirable function. RNA-based reg- 
ulation is an alternative approach for controlling gene 
expressions [T5j [Sj [Mj |34] . RNA structures and inter- 
actions with other nucleic acid species can be reliably 
predicted based on Watson-Crick base-pairing, much 
more so than typical protein-protein or protein-DNA 
interactions of protein regulators. As for the choice of 
biochemical infrastructure, the unintended interactions 
between the circuit and its environment can be greatly 
reduced by reconstructing the circuit in vitro. In vitro 
implementation of efficient transcription and transla- 
tion machinery |30 [ I17 ) for synthetic networks have been 
successfully implemented |26] . Yet, a supporting envi- 
ronment for an in vitro RNA-based regulatory circuit 
can be even simpler as there is no need for translation, 
protein maturation, and protein-DNA interactions. 

Previous work |191 120] introduced in vitro transcrip- 
tional circuits as simplified synthetic genetic regulatory 
circuits. Individual switches functioning as inverters 
and a bistable feedback circuit composed of two in- 
verters have been demonstrated (Figure lA, top). Our 
"DNA switch" , a simplified gene, has a promoter for T7 
RNA polymerase (RNAP) fianked by two separate do- 
mains, an input domain and an output domain. Down- 
stream of a promoter lies an output domain that en- 
codes "an RNA product" ; on the opposite side of the 
promoter, an input domain regulated by "an RNA reg- 
ulator" via simple Watson-Crick base-pairing rule is lo- 
cated. The modularity of a DNA switch allows for an 



independent design of an RNA product and an RNA 
regulator within a switch. Hence, one can "wire" sev- 
eral switches together to compose a complex regulatory 
network, in principle, by simply designing the RNA 
output of one switch to be the RNA regulator of the 
other switch. Moreover, individual switch characteris- 
tics such as switching thresholds and maximum output 
levels are set by the concentrations of switch compo- 
nents rather than by molecular characteristics of bind- 
ing domains. The state of each switch (transcription 
rate) and the levels of signals (RNA concentrations) re- 
layed among switches define the behavior of the overall 
circuit. As the circuit dynamics relies only on RNA 
transcription and degradation, our in vitro circuit op- 
crates in a relatively simple environment with NTP fuel 
and only a few enzymes, RNA polymerase and RNases. 

Here, we expand the repertoire of circuit motifs for 
transcriptional circuits by implementing a repeater. A 
repeater is a transcriptional switch whose RNA output 
level is a sigmoidal activation function of its RNA reg- 
ulator, thus providing a concise mechanism for relaying 
activation signal. An activator offers greater fiexibil- 
ity and simplicity for a synthetic circuit design and al- 
lows for faster timing. Although, computationally, two 
inverters connected in a series can substitute for a re- 
peater, such design can lead to delay and unintended 
amplification of noise [13] . In addition, a repeater wired 
to itself can implement a positive autoregulatory switch 
that simplifies our previous bistable circuit design |20] 
as demonstrated in synthetic in vivo circuits [H 116] 
(Figure lA, bottom). In this report, we describe the 
design for a repeater switch and characterize each ele- 
mentary reaction required for a functional repeater. We 
then construct and analyze a self-activating switch, a 
single repeater wired to itself, that can be tuned to ex- 
hibit bistability. Furthermore, we characterize its sen- 
sitivity to parameters such as DNA and enzyme con- 
centrations; we show that the experimental results are 
consistent with a simple mathematical model, whose 
parameters are set by Bayesian inference from a subset 
of experimental data. 

Inverter versus Repeater 

The repeater that we implement in this paper shares 
several components and functional mechanisms with 
the previously implemented inverter. Therefore, we 
first describe the components and functional mecha- 
nisms of a transcriptional inverter. An inverter con- 
sists of two components, a DNA template ("T") and 
a DNA activator ("A"). The DNA template T consists 
of single-stranded regulatory domain, partially single- 
stranded T7 RNAP promoter, and a double-stranded 
region encoding an RNA output ("rO"). This partially 
single-stranded promoter is transcribed poorly by T7 
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RNA polymerase [23l [20|, and thus, is designated as 
an OFF state. The single-stranded DNA activator A 
is complementary to the missing promoter region of T. 
Upon hybridization of T and A, the resulting T-A com- 
plex has a complete promoter except for a nick and was 
found to be transcribed well, approximately half as effi- 
ciently as a fully double-stranded promoter. Therefore, 
T-A is designated as an ON state. The inverter is reg- 
ulated by an RNA inhibitor ("rl") that is complemen- 



tary to A. Because A-rl complex is thermodynamically 
more stable than T-A complex and the single-stranded 
domain of A beyond the helical domain of T-A complex 
is available for initiating hybridization reaction with rl, 
rl can strip off A from T-A complex through a toehold- 
mediated branch migration reaction [25l [28] . Typically, 
A is in excess of T such that the input rl will first re- 
act with free A, then strip off A from ON-state switch, 
and the remaining rl will be free-fioating in solution. 
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Figure 1: Transcriptional switches and circuits. (A) Schematic diagrams of an inverter, a repeater, and circuits. 
Blunt ends indicate inhibition, while arrowheads indicate signal production or activation. (B) An inverter switch. T 
is the DNA template with an incomplete promoter region, A is the DNA activator, rl is the RNA inhibitor, and rO 
is the RNA output. Sequence domains are color-coded such that the same colors represent either complementary or 
identical sequences. Hybridization reactions are marked by black arrows. Transcription by RNAP and degradation by 
RNase H and RNase R are marked by black dashed arrows and dotted arrows, respectively. The production curve of 
rO as a function of rl, the degradation curve of rO as a function of rO, and the steady-state curve of rO as a function 
of rl constructed from the composition of the production and degradation curves are shown. Green dashed lines mark 
[rl'°*] — [A'°'] - [T'°*], where all free A species are consumed by rl, yet the switch is still fully ON; orange dashed 
lines mark [rl'°'] — [A*°'], where all A species are consumed such that the switch is fully OFF. The purple dashed line 
marks [rO'°'] — [dX*°'] where dX is the regulatory target of rO: below this level rO is mostly bound to dX such that 
degradation curve is dominated by RNase H, while above this level RNase R degrades free rO. A dashed red line in 
degradation plot illustrates the case when only RNase H is present in the reaction mixture. Note that the annihilation 
reaction between A and rl is not shown (see Figure 3). Also, dX strand or RNase reactions on rO are not shown. (C) 
A repeater switch. T and A are the same as described for an inverter switch, however, a repeater switch also uses a 
DNA inhibitor dl to set the switching threshold. Green dashed lines mark [rA*°'] = [dl'°*] - [A*°'], where all free dl 
species are consumed by rA, yet the switch is fully OFF; orange dashed lines mark [rA*°'] — [dl'°*] - [A'°*] + [T*°'], 
where enough A species are freed from dl and available to hybridize with T, such that the switch is fully ON. The 
purple dashed line marks [rO'°'] = [dX'°*], where dX is the regulatory target of rO as above. Note that annihilation, 
interfering, recovering, recapturing reactions, dX strand, and RNase reactions on rO are not shown (see Figure 3). 



Overall, the production rate, as well as the fraction of 
ON-state switch, is a sigmoidal inhibitory function of 
rl (Figure IB, left). 

The degradation speed of RNA signals plays an im- 
portant role in circuit dynamics by setting the time- 
constants of signal relays. At the same time, the shape 
of the degradation curve in combination with the pro- 
duction curve would determine the steady-state RNA 
levels. The degradation function is determined by the 
concentrations of substrates and the enzyme constants 
of RNases, Escherichia coli ribonuclease H (RNase H) 
and RNase R. RNase H degrades RNA that is hybridized 
to DNA: rl in A-rl complex and rO that is bound to its 
downstream regulatory target, "dX", to form a dX-rO 
complex. At low rO concentrations, most of rO species 
exist within dX-rO complex such that the degradation 
rate is largely dictated by RNase H. On the other hand, 
RNase R degrades single-stranded RNA, both free rl 
and free rO, such that the degradation rate of free 
rO — the amount of rO in excess of total dX con- 
centrations — will be largely determined by RNase R. 
Because RNase H has a low Michaelis constant, the 
degradation rate of rO by RNase H quickly saturates as 
the concentration of dX-rO increases. Also, the degra- 
dation rate of rO by RNase R saturates as free rO 
concentration increases. Each degradation curve by 
RNase H and RNase R constitutes a typical Michaelis- 
Menten saturation curve with different origins: [rO*°*] 
= for RNase H and [rO*°*] = [dX'°*] for RNase R. 
Thus, the composition of degradation curves by these 
two enzymes results in the degradation function with 
a kink located at the total concentration of dX (Fig- 



ure IB, middle). When the maximum rates and switch- 
ing thresholds for production and degradation curves 
are approximately matched, the resulting steady-state 
RNA output level shows a sigmoidal inhibitory response 
with respect to RNA inputs, although the transition re- 
gion contains a kink (Figure IB, right). 

A repeater utilizes some of the same modular design 
motifs as an inverter. Because RNA polymerase tran- 
scribes poorly from a DNA/RNA hybrid promoter 
we chose to implement a repeater through an indirect 
activation by an RNA activator. The transcriptional 
repeater also contains a DNA template T and a DNA 
activator A. However, the repeater has an extra com- 
ponent, a DNA inhibitor ("dT'). Much like the RNA 
inhibitor rl of an inverter, the DNA inhibitor dl can 
bind to and remove A from the ON state T-A complex 
and form a A-dl complex. The input of the repeater, 
an RNA activator. "rA" . is a single-stranded RNA that 
can displace dl within the inhibiting complex A-dl and 
release A through a toehold-mediated branch migra- 
tion. Then, the released A can bind back to T and turn 
the switch on. Unlike the inverter, the concentrations 
of T and A are about the same and the concentration 
of dl is in excess of A to provide activation thresh- 
old for rA. (The DNA activator strand concentration 
should be roughly comparable with the template con- 
centration; it should be at least as high, so that all the 
template can be turned ON, but it needn't be higher, 
since excess activator merely disables a stoichiometric 
amount of inhibitor.) The input rA will first react with 
free dl, then strip off dl from the A-dl complex, and the 
remaining rA will be free-floating in solution. Overall, 
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the production rate, as well as the fraction of ON-state 
switch, is a signioidal activation function of rA (Fig- 
ure IC, left). The degradation function of a repeater 
is identical to that of an inverter (Figure IC, middle). 
Therefore, the resulting steady-state RNA output level 
shows a sigmoidal activation response with respect to 
RNA input with a kink within the transition region 
(Figure IC, right). 

Results 

System Design and Mechanisms 

The challenge of designing sequences for a functional 
repeater lies in accommodating desirable hybridization 
reactions outlined above and suppressing side reactions 
and crosstalks. Here, we present the sequence design 
for a repeater and provide evidence for its proper func- 
tionality. 

Sequences of the repeater components are chosen to 
minimize alternative folding [S] and spurious interac- 
tions [29]. Utilizing the modularity of synthetic switch 
designs, several domain lengths have been adapted from 
previous work |20| : for example, the binding domains 
of an OFF switch template to an activator (27 bases) 
and the toehold of an activator (9 bases). The 3' end 
hairpin structure (16 bases) of RNA output increases 
copy number and also decreases self-coded extension of 
RNA transcripts by RNAP [33]. The sequences of T, A, 
dl, rA, and their complexes are shown (Figure 2). Note 
that, in this design, the sequence of RNA output, rO, 
is identical to RNA input, rA, because we constructed 
and characterized a self-activating switch. Because we 
implemented an indirect activation by RNA inputs, the 
RNA activator rA and DNA activator A shares common 
sequence domains. Thus, a T-rA complex is expected 
to form when both T and rA are available. However, 
it is not desirable if an excess of rA interferes with the 
hybridization reaction between T and A. Thus, we im- 
plemented a staggered design, where the A-dl complex 
leaves 5 bases of A as single-stranded (colored light blue 
in Figure 2) and the rA-dl complex leaves 4 bases of 
dl as single-stranded (colored dark blue in Figure 2), 
resulting in a total of 9 base-pair differences between 
the proper ON state, T-A, and the interfering complex, 
T-rA. We chose the sequences and lengths of the bind- 
ing domains for T, A, dl, and rA such that the pre- 
dicted AG° of complexes are in the following order: 
rA-dl < A-dl < T-A < T-rA. In terms of the number of 
base pairs, there are 18, 27, 31, and 34 complementary 
base-pairs for T-rA, T-A, A-dl, and rA-dl complexes, 
respectively. Therefore, the RNA activator rA has the 
strongest influence on the state of a repeater by design. 
This design meets the requirement of proper hybridiza- 



tion reactions for a repeater as discussed below. 

There are four types of simple hybridization reac- 
tions which we call activation, annihilation, subduing, 
and interfering (Figure 3A, red boxes): A binds to T, 
forming T-A (activation); dl binds to A, forming A-dl 
(annihilation); rA binds to dl, forming rA-dl (annihila- 
tion); rA binds to T, forming T-rA (interfering). Simple 
hybridization reactions are thermodynamically favor- 
able because the resulting complex gains several base 
pairs, and hence the reactions proceed in a unidirec- 
tional way. As mentioned above, 'interfering' reaction 
is not desirable and we provide kinetic pathways (re- 
covering and recapturing) to quickly resolve the inter- 
fering complex as shown below. There are four types of 
strand displacement reactions, which we call inhibition, 
release, recovering, and recapturing (Figure 3A, orange 
boxes): dl strips off A from T-A complex to form A-dl 
(inhibition); rA displaces A from A-dl complex to form 
rA-dl (release); A displaces rA from T-rA complex to 
form T-A (recovering); dl strips off rA from T-rA com- 
plex to form rA-dl (recapturing). All displacement re- 
actions are designed to be initiated at the 'toehold,' a 
single-stranded overhang that extends beyond the heli- 
cal domain of the initial complex. The incoming strand 
can bind to this toehold, providing a fast kinetic path- 
way through a toehold-mediated strand displacement 
reaction [551 HH] • The predicted thermodynamic ener- 
gies of starting complex and resulting complex indicate 
that all displacement reactions are thermodynamically 
favorable and hence approximately unidirectional. 

The hybridization reactions and strand displacement 
reactions for switch components are characterized by 
running different combinations of T, A, dl, and rA 
in a non-denaturing gel (Figure 3A). The two single- 
stranded species comprising T was annealed prior to 
mixing with other single-stranded species. The indi- 
cated components were simultaneously mixed in a test 
tube at room temperature and were allowed to sit for 
five minutes before being subjected to a non-denaturing 
gel run at 4°C. When T and A were mixed together, A 
could bind to T, resulting in a well-defined band of T-A 
complex (lanes 1 vs. 2: activation). When T, A, and 
an excess of dl were mixed together, only T and A-dl 
complex were observed, implying that dl can bind to A 
and also strip off A from T-A (lanes 2 vs. 3: annihila- 
tion and inhibition). Of note, the single-stranded A and 
dl migrated faster than the 30-bp marker and were not 
visualized in the gel. Also, when rA and T are mixed to- 
gether, rA can bind to T and form a interfering complex 
T-rA, which migrated slower than T-A complex (lane 
4: interfering). However, when all four components 
were mixed together with the sum of A and rA concen- 
trations in excess of dl, T-rA complex was no longer 
visible, while T-A and rA-dl complexes were observed. 
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T 

CATTAGTGTCGTTCGTTCACAGTAATA CGACTCACTATAGGGAGAAACAAAGjVACGAACGACACTT^TGAACTACTACTACACACTjVATACTGACAT^GTCAGAAay. 




T-A 

CATTAGTGTCGTTCGTTCACAGTAATA CGACTCACTATAGGGAGAAACAAAGAACGAACGACACTAATGAACTACTACTACACACTAATACTGACAAAGTCAGAAAy 
\DVID¥IDWDJ,WJ,DVD¥DDWDDW0IDJ,DVJ.lVJ,\3DJ-3¥DJ,0¥J,¥J,JDDJ,JiJ,XDJ,I,lDI,13DlJ,DDJ,D10¥lJ,¥DJ,J,DVJ,DVJ,0VJ,D10I,3¥I,I¥I,0¥DlDJ,XJ,D¥DJ,DiJ.i 




Figure 2: Sequences of DNA and RNA species of a self-activating switch. The sequence domains are color-coded to 
indicate identical or complementary sequences. Although we simply colored the output domain of switch coding for 
regulatory output ('rO' in Figure 1) red, this sequence is identical to the RNA input 'rA' because we constructed a 
self-activating switch. 



implying that A and dl can strip off rA from T-rA and 
dl can bind to rA (lanes 4 vs. 5: recovering, recaptur- 
ing, and annihilation). Moreover, it is also implied that 
rA can bind to dl and strip off dl from the A-dl complex 
and release A, which in turn binds to T, resulting in a 
T-A complex (lanes 3 vs. 5: annihilation, release, and 
activation). We did not observe a noticeable A-dl band 
in lane 5, indicating that the concentration of rA may 
have been underestimated for lane 5. Recent reported 
measurements of DNA hybridization rates ranged from 
10^ to 10^/M/s [10]. Thus, starting from equal concen- 
trations of switch components at hundreds of nM, the 
hybridization reactions would be completed with half- 
lives on the order of 10 to 100 seconds. No smearing 
or reaction intermediates were detected between the gel 
bands when several switch components were mixed to- 
gether, indicating that the hybridization reactions and 
strand displacement reactions were mostly completed 
within a few minutes. 

Three types of enzyme reactions are separately char- 



acterized as follows (Figure 3B). RNAP could efficiently 
transcribe rA from an ON-state template, T-A. How- 
ever, the transcription was much slower from an OFF- 
state template, T (Figure 3B, top). Thus, transcription 
reaction is much more efficient from an ON-state switch 
than from an OFF-state switch as desired. In this case, 
since dl was not provided in the transcription reaction, 
we do not expect to observe auto-regulation. RNase H 
could degrade rA when both rA and dl were present, 
but no degradation of rA was observed without dl (Fig- 
ure 3B, middle). Thus, RNase H degrades RNA that 
are hybridized with DNA, but not free-floating RNA. 
RNase R could degrade about 4 fiM rA within 210 min- 
utes. On the other hand, when 4 of rA was mixed 
with 1 fiM of dl, about 1 /iM of rA was left over af- 
ter 210 minutes (figure 3B, bottom). Thus, RNase R 
degrades single-stranded RNA, but not RNA within 
RNA-DNA hybrid complexes. 

Taken together, our sequence design led to proper 
hybridization and strand displacement reactions among 
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switch components with fast kinetics. The enzyme re- 
actions could provide production and degradation of 
RNA regulatory signals with specific recognition of sub- 
strates. Also, note that the time-scales of enzyme re- 
actions are much slower than typical hybridization or 
strand displacement reactions. Therefore, the presumed 
sharp thresholds achieved by fast and irreversible hy- 



bridization kinetics would be approximately valid, even 
in the presence of constant enzyme-mediated produc- 
tion and degradation of signals. 

Self-activating switch 

To demonstrate the functionality of the repeater 
mechanism, we chose to implement — with a single 
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Figure 3: Schematic representation of reactions for a repeater switch. (A) DNA and RNA hybridization reactions. The 
red and orange boxed reactions show simple hybridization reactions and strand displacement reactions, respectively. 
Note that, all reactions are thermodynamically driven and that toehold-mediated branch migrations provide fast 
kinetic pathways for the strand displacement reactions. A non-denaturing gel was used to analyze the results of such 
hybridization and strand displacement reactions in the absence of enzymes (top right). The DNA and RNA species 
were mixed and incubated for 5 min at room temperature prior to loading onto the non-denaturing gel. For each lane, 
the red label indicates the constituents loaded in that lane: T is 50 nM of [T'°*]; A is 500 nM of [A'°*]; dl is 700 
nM of [dl*°']; and rA is 500 nM of [rA*°']. The leftmost lane contains a 10 base-pair ladder. The blue labels on the 
right side of gel indicate corresponding single-stranded species and complexes. Single-stranded A and dl bands do not 
appear because they have been run off of the gel. (B) Enzyme reactions. The reaction diagrams are shown in blue 
boxes with the corresponding denaturing gel analysis of experimental results on their right sides. For the denaturing 
gels, the leftmost lanes contain 10 base ladders and the templates or substrates for enzymes are indicated by red labels 
on top of gels. Samples in lanes 1 through 4 (and also for lanes 5 through 8) were taken from the same reaction tube 
at different time points (0, 10, 40, and 210 min). Sample in lane 9 of the gel analyzing RNAP reaction contained a 
half volume of that loaded in lane 8 to avoid SYBR gold signal saturation. Control lanes (lanes 10 and 11 for the gel 
analyzing RNAP reaction and lanes 9 through 11 for gels analyzing RNase H and RNase R reactions) have purified 
rA at the indicated concentrations. Enzyme concentrations used are as follows; [RNAP] = 33.4 nM, [RNaseH] = 
0.168 nM, and [RNaseR] = 0.336 nM. 



self-rcgulating repeater switch — a bistable latch whose 
function is similar to that of the two-node network de- 
veloped in [5D]. The construction of a self-activating 
switch is straightforward as it only requires the RNA 
output rO to be identical to the RNA input rA, as 
shown in Figure 2. Although we have shown that the 
elementary reaction steps are plausible, the repeater 
function embedded in a transcriptional circuit may show 
dynamic features deviating from our expectation. For 
instance, previous works demonstrated that the switch 
behavior was measurably different when driven by an 
RNA species subject to constant production and degra- 
dation as compared to being driven by a DNA species [5D] 
Therefore, the motivation is to assess how well the self- 
activating switch behaves with respect to its expected 
behavior and to gain quantitative understanding of ac- 
tivation switch motif within a system. Further, this 
simple circuit highlights the simplicity provided by the 
activation switch motif vis-a-vis creating the same func- 
tionality from a larger network of inhibitory switches. 
In the following sections, we first 'qualitatively' ex- 
plain the dynamics of self-activating switch based on 
the production and degradation functions of RNA sig- 
nals. Then, we assess how well the kinetic model using 
ensemble parameters matches the experimental results 
quantitatively. 

Self-activating switch: qualitative and quantita- 
tive prediction of bistability 

For a self-activating switch, the RNA input rA and 
RNA output rO are identical, and therefore, we call 
RNA output as 'rA' henceforth. The target DNA of the 
RNA output, dX, is now the DNA inhibitor, dl. There- 



fore, the production and degradation curves shown in 
Figure IC can now be interpreted as the production and 
degradation curves of RNA activator rA as functions of 
its own concentrations. Thus, given the same switch 
template and DNA activator concentrations ([T'°*] = 
[A*°']), the concentration of rA that consumes all dl 
and turn the switch ON in the production curve coin- 
cides with the threshold in the degradation curve below 
which degradation by RNase H is dominant (i.e. [rA*°'] 
= [dl*°']). The intersections of production and degra- 
dation curves thus indicate the steady-states where to- 
tal rA concentrations remain constant, and where and 
how they cross will dictate whether the self-activating 
switch show monostable or bistable behavior. If an un- 
stable fixed point exists where a slight increase in rA 
drives further rA production and a slight decrease in 
rA drives further rA reduction (Figure 4 A, right), the 
system will show bistable behavior with the threshold 
set at this unstable fixed point. In that case, the switch 
states approach either of the two stable steady-states, 
completely ON or OFF, depending on the initial con- 
centration of rA. On the other hand, a monostable ON 
state will be achieved irrespective of initial RNA acti- 
vator concentration if the production rate exceeds the 
degradation rate at the switching threshold (Figure 4A, 
left); a monostable OFF state will be achieved irrespec- 
tive of initial RNA activator concentration if the degra- 
dation rate exceeds the production rate of the RNA ac- 
tivator when the switch is fully ON (Figure 4A, center). 
The bistable and monostable switch behaviors depend 
not so much on the physical nature of molecular com- 
ponents but on the features set by continuously tunable 
concentrations — such as the output amplitude and the 
activation threshold. 
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Figure 4: Bistability of a self-activating switch. (A) Production and degradation functions, and dynamics of a self- 
activating switch. Three plots illustrate possible configurations of production and degradation functions of rA as a 
function of rA. Depending on the locations and stabilities of fixed points, two types of monostable behaviors or a 
bistable switch behavior are expected. (B) A denaturing gel analysis to test bistability of a self-activating switch. 
Lanes 1 through 12 show the results from 12 separate reactions. The reaction conditions are as follows: [T*°*] = 50 
nM, [A*°'] = 90 nM, [dl'°*] = 1000 nM, [RNAP] = 66 nM, [RNaseH] = 0.7 nM, and [RNaseR] = 0.23 nM, with 

variable initial rA concentrations ([rA'°*] — 0, 0.1, 0.2 1, and 3 fj,M). After 120 min, all reactions were stopped and 

subject to denaturing gel analysis. Control lanes, lanes 13 and 14, have purified rA at the indicated concentrations. 
Brackets marked 'w' indicate incomplete transcription products and degradation products. (C) The measured final 
concentrations of rA from the denaturing gel in (B). Note that the apparent threshold of switching between the low 
and high steady-states is between 0.6 and 0.7 /iM of initial rA concentrations. Black dashed line marks the points 
where the initial and final rA concentrations coincide. 



We then experimentally tested whether the self-acti- 
vating switch can show bistable output responses (Fig- 
ure 4B). When the self- activating switch was subject 
to transcription and degradation reactions with a wide 
range of initial RNA activator concentrations, the fi- 
nal RNA activator concentrations reached two distinct 
states after two hours: almost nM or higher than 
600 nM. The self-activating switch showed bistable be- 
havior as expected. However, the experimental results 
deviated from an idealized circuit behavior in two ways: 
a low switching threshold and non-constant ON state 
outputs. First, the switching threshold lied between 600 
and 700 nM of rA in the experiment while the expected 
threshold from idealized production and degradation 
curves lied between 910 and 960 nM of rA for the case 
shown in Figure 4B. This discrepancy may be explained 
by the "burst phase" in enzyme kinetics [TH], which 



could add a few copy-numbers of RNA transcripts rA 
and effectively lowered the apparent switching thresh- 
old for initial rA concentrations. Second, the ON state 
output levels were not constant in contrast to almost 
identical OFF state output levels. High final RNA 
outputs were measured for high initial RNA concen- 
trations, which could be explained by a slow degrada- 
tion kinetics near the steady-state. On the other hand, 
the measured RNA output levels were as low as 600 
nM when more than 960 nM of rA would be required 
to maintain a fully ON state in the idealized model. 
We used RNase R to clean up the incomplete single- 
stranded degradation products by RNase H. Neverthe- 
less, incomplete degradation products accumulated over 
time (Figure 4B, brackets marked 'w'). Thus, it is pos- 
sible that some portions of these incomplete degrada- 
tion products still function much like a full-length RNA 
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activator — e.g. rA without the hairpin region at the 
3' end. 

Due to experimental difficuhics including RNase R 
not being commercially available at the time of experi- 
ment and having apparently short life-time, wc focused 
our study on a simpler system without RNase R. Fur- 
ther, excluding RNase R eliminates the kink within 
degradation function and possibly allow for a simpler 
quantitative interpretation. However, since the degra- 
dation curve by RNase H is insensitive to the change 
of rA concentrations when the concentration of rA is 
above the total concentration of dl, the production and 
degradation curves may not cross at the high steady- 
state. In that scenario, the high state of rA concen- 
tration can grow without bound — which make the 
system behavior more latch-like. Even so, the behavior 
of switch state is bistable. 

To test bistability of switch response in the absence 
of RNase R, we monitored the kinetics of switch re- 
sponse through real-time fluorescence measurement. The 
3' end of A is labeled with a Cy5 dye and the 5' end 
of dl is labeled with an lowaBlack-RQ quencher such 
that the fluorescence is low when the A-dl complex is 
formed due to fluorescence quenching [22] while the flu- 
orescence is high when A is free or within the T-A com- 
plex (Figure 5A). Of note, fluorophore-quencher inter- 
action can stabilize the resulting complex, A-dl 
(In preliminary work, we observed mild to severe dif- 
ferences between fluorophore-labeled strands and their 
unlabeled twins, depending on the fluorophore used.) 
We initiated the reaction with four different rA con- 
centrations and measured how the system approached 
different steady-states (Figure 5B). The fluorescence 
traces either reached a maximum signal, which implied 
a completely ON state, or a minimum signal, which im- 
plied a completely OFF state. The fluorescence moni- 
toring indicated that the self-activating switch quickly 
approached steady-states and that both fully ON and 
OFF cases were stable steady-states. At the same time, 
the dynamics of RNA signal was determined by taking 
samples at different time points and measuring the to- 
tal rA concentrations in gel (Figure 5CD). The high 
fluorescence traces (magenta and blue) corresponded 
to more than 2 /iM final RNA levels, but the low fiuo- 
rescence traces (green and red) corresponded to about 
0.5 /iM final RNA levels. Thus, we achieved bistable 
behavior in both switch and RNA signals, albeit the 
fact that the high RNA output level is not bounded. 

While the bistable system behavior was demonstrated 
experimentally, a mathematical framework for the in 
vitro self-activation system would be required to inves- 
tigate several further questions. First, do the elemen- 
tary reactions prescribed in Figure 3 explain the system 
behavior? Second, how fine-tuned the rate constants 



have to be to achieve bistability? Third, are the esti- 
mates and bounds for rates obtained from the char- 
acterization of elementary reactions compatible with 
the system behavior? To address these questions, we 
constructed a simple mathematical model for the self- 
activating switch that uses four ordinary differential 
equations. The dynamics of each DNA and RNA species 
were derived from the hybridization and strand dis- 
placement reactions described in Figure 3, assuming 
Michaelis-Menten enzyme kinetics (see Appendix) . Note 
that the concentration of 'interfering complex' T-rA was 
assumed to be negligible in the simple model, and thus, 
'interfering', 'recovering', and 'recapturing' reactions 
were not included. Initial simulation results showed 
that a bistable switch response can be achieved with 
plausible rate parameters similar to a previous work [20] 
(data not shown). 

To test whether our mathematical model can be eas- 
ily generalized using a relatively small number of exper- 
iments to determine system behavior for a wider range 
of parameter choices, we need a separate training and 
test data sets. Yet, a rich enough training data set 
would be required to constrain parameters effectively. 
Thus, the training data set was obtained by measur- 
ing final rA concentrations for self-activating switches 
initialized with a wide range of rA concentrations and 
three different dl concentrations as switching thresholds 
(Figure 6A, red, blue, and green). For this data set, fiu- 
orescence data for available A concentrations were not 
measured. Below the apparent switching threshold, the 
final concentrations of rA approached ^0.5 /iM regard- 
less of the initial rA concentrations. On other other 
hand, for initial rA concentrations higher than the ap- 
parent switching thresholds, the final rA concentrations 
were distinctively higher. As expected from the analysis 
of idealized production and degradation curves shown 
in Figure 5A, adjusting dl concentrations shifted the 
threshold of switching accordingly. 

To find suitable values for the set of 11 parameters 
in our model, we used a Monte Carlo Bayesian inference 
approach that results in an ensemble of parameter sets 
that are compatible with the given data within noise 
bounds [S]. (See Appendix for details.) After fitting 
to the experimental results shown in Figure 6A, the re- 
sulting ensemble of parameter sets can be used to make 
ensemble predictions for novel experimental conditions; 
for example, we can estimate the possible range of con- 
centrations of any chemical species at any time point 
given initial conditions. 

Using the ensemble parameters fitted to the train- 
ing set (Figure 6A), the simulation results for Figure 5 
were generated and compared with the experimental 
results (Figure 6B). The kinetic model predicted that 
the switch would either approach the low steady-state 
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Figure 5: Time courses of a self-activating switch. (A) The schematic reaction diagram of a self-activating switch 
(cf. Figure IC). The activator A is labeled with Cy5 fluorophore (red circle) and the inhibitor dl is labeled with 
lowaBlack-RQ quencher (black circle) for real-time monitoring of switch states by fluorescence. A possible arrangement 
of production and degradation curves for self-activating switch without RNase R in the reaction mixture is shown as an 
inset (cf. Figure 4A). (B) For fluorescence measurements, all the switch components, DNA template T, fluorophore- 
labeled A, and quencher-labeled dl, were added to each cuvette. Then, purified rA was added prior to enzyme 
additions such that the initial conditions are different for each cuvette. Finally, the enzymes were added at time 
and thoroughly mixed (black dashed line). Normalized fluorescence signal measures the fraction of A that is not 
bound to dl, i.e. ([A] + [T A])/[A'°*], because the fluorescence from A is quenched when A is in a A-dl complex. All 
four cuvettes contained 48 nM of [T*°'], 145 nM of [A'°*], 1500 nM of [dl*°'], 16.7 nM of [RNAP], and 1.68 nM of 
[RNaseH] with the initial rA concentrations at 0, 550, 1350, and 1750 nM (colored red, green, blue, and magenta, 
respectively). (C) Time courses of [rA*°'] with samples taken from four different cuvettes. The colors correspond to 
those of fluorescence traces in (B). (D) One of the denaturing gels for [rA'°'] measurements as shown in (C). Samples 
in lanes 1 through 5 (as well as lanes 6 through 10) were taken from the same cuvette at different time points. The 
final rA concentrations were measured with respect to the known concentrations of purified rA in the control lanes, 
lanes 11 through 13. 



in which ah A were bound to dl and the rA concen- 
trations converged to 0.5 or the high steady-state 
in which ah A were free from dl and the rA concen- 
trations kept increasing. The experimentally measured 
fluorescence signal from A and sampled rA concentra- 
tions were close to the prediction (Figure 6B), despite 
the fact that the training data set contained informa- 
tion on rA concentrations at a single time point only. 
However, in the experiment, the rA concentrations for 
the ON state did not increase indefinitely as the model 
predicted, but seemed to reach a high steady-state. We 
shall address this discrepancy between the model and 
experiment later. 



To further test the sensitivity of system behavior to 
parameter variation, we systematically varied the tem- 
plate ([T'°*]), activator ([A*°*]), inhibitor ([dl*°']), and 
enzyme concentrations ([RNAP] and [RNaseH]) (Fig- 
ure 6C). For each parameter choices, we initiated the 
transcription and degradation reaction with either a 
low amount of rA (0 nM) or a high amount of rA 
(730 nM). If the final difference of RNA activators is 
greater than the initial difference of 730 nM, it would 
imply that these two initial conditions lead to differ- 
ent switch states. However, since the switch states 
were not monitored by fluorescence, these results were 
not a direct readout of bistable switch response. On 
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Figure 6: Quantitative measurements and model prediction of the self-activating switch dynamics. The circles within 
each figure and the lines with darker colors in (B, left) are from experimental measurements. Lines with lighter 
shades are generated from 95% confidence interval of model parameters fitted to the experimental results shown in 
(A). Experimental results in (B) and (C) are used as test cases of model prediction. (A) The training data set used 
for parameter fitting. The circles are gel measurements of final [rA'°*] at 210 min. All reactions contained 48 nM 
of [T*°'], 145 nM of [A'°*], 16.7 nM of [RNAP], and 1.68 nM of [RNaseH]. The concentration of dl was varied as 
follows: [dl'°*] = 1.13 /iM (red), 1.45 /iM (blue), or 1.87 ^M (green). Note that the thresholds of switching between 
low and high states were apparently shifted as the total concentration of dl was increased. (B) Experimental and 
predicted time courses of switch states (monitored as fluorescence changes) and rA concentrations. The experimental 
data is identical to Figure 5. The experimental conditions were similar to those colored blue in (A): only dl was 
50 nM higher for (B). (C) Experimental and predicted parameter dependence of self-activation switch behavior. The 
default set of experimental conditions were identical to those colored red in (A): [T*°'] = 48 nM, [A'°*] = 145 nM, 
[dl*°'] = 1.13 ^iM, [RNAP] = 16.7 nM, and [RNaseH] = 1.68 nM (marked as black dashed lines in each subplot). 
For each experimental measurement, one parameter was varied at a time and the final concentration of rA after 210 
min was recorded with initial rA concentration of nM (red) or 730 nM (blue), respectively. 



the other hand, if both initial conditions lead to low 
steady-states, the switch could be in the monostable 
regime, but it also could be in a bistable regime that 
requires higher initial rA inputs to reach high steady- 
state. Yet another possibility is when both initial condi- 
tions lead to high steady-states such that the differences 
of final rA concentrations would be close to 730 nM. 
Thus, for experimental results, we marked out the pa- 
rameter ranges where the differences of final rA con- 
centrations were greater than 800 nM as indications of 
bistable switch responses. The model predicted that 
the self-activating switch showed a bistable response 
for a limited range of parameter variations near the 
default values. For simulation results, we map the pa- 



rameter range as bistable, if the final switch states were 
different for the 95% of ensemble parameter values. 
From the simulation results, the self- activating switch 
reached the OFF state under conditions with low T, 
low A, low RNAP, high dl, or high RNase H, irrespec- 
tive of initial conditions. Conversely, the self-activating 
switch reached the ON state under conditions with high 
T, high A, high RNAP, low dl, or low RNase H, irre- 
spective of initial conditions. The experimental results 
were in good agreement with the model prediction on 
the range of bistability for each parameter. Not sur- 
prisingly, the model predictions showed increased un- 
certainty in the transition regions for both high and 
low initial rA concentrations. 
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Discrepancies in modeling 

Given its simplicity, our model predicted the dy- 
namics of the switch surprisingly well. We only in- 
cluded the main hybridization and enzyme reactions 
that constituted our design intentions, excluding the 
unintentional 'interfering,' 'recovering,' and 'recaptur- 
ing' reactions as well as any other undesired crosstalk 
between DNA or RNA species. Notably, the training 
data used for fitting (Figure 6A) only gives informa- 
tion about how the final [rA'°*] depends on the initial 
[rA*°*] and [dl*°'], and thus provides very weak kinetic 
constraints on the model. Unsurprisingly, a wide range 
of parameter sets could fit the data within reasonable 
noise bounds, as is common when fitting biochemical 
system models to limited data [TT]. Nonetheless, qual- 
itative (and sometimes quantitative) predictions were 
made accurately for a variety of novel experimental con- 
ditions (the test data. Figure 6B,C), including the time 
evolution of the RNA signal and switch state. This 
suggests that the intended reaction pathways, used to 
construct the model, are dominant ones in the system 
dynamics. 

However, some predictions were qualitatively wrong 
for all parameter sets in the a posteriori ensemble. These 
prediction errors may point to omissions in the model. 
As an example, for the time evolution of the RNA signal 
(Figure 6B, right), the model predicted that the upper 
two reactions would reach a state in which [rA'°*] keeps 
increasing, while in the experiment, these two reactions 
seemed to converge to a steady-state. A possible expla- 
nation is that incomplete degradation by RNase H, as 
observed in previous work |20] . yields short fragments 
of rA that accumulate during the reaction: at high con- 
centrations of these short products, product rebinding 
to the transcription initiation complex of RNAP may 
decrease transcription rates [21]. 

As another example, the model predicts a stronger 
sensitivity to RNase H concentration than was observed 
(Figure 6C, rightmost panel). Incomplete degradation 
fragments may also provide an explanation for this dis- 
crepancy. These shorter oligonucleotides could also bind 
to dl and thus compete for the degradation by RNase H, 
decreasing its effectiveness. Other notable discrepan- 
cies are that the model predicts that wider ranges of 
[T*°*] and [A*°'] wiU sustain bistability than was exper- 
imentally observed. Another oversimplification of the 
model is that the inhibition and release reactions are 
expected to be reversible, despite being thermodynam- 
ically favored in the forward direction as described; they 
are instances of so-called toehold exchange, as there 
is also a (weaker) toehold in the reverse direction as 
weU [36]. 



Conclusions 

In this work we have designed and synthesized an 
activating cascade for in vitro transcriptional circuits 
and demonstrated its use as a single-switch bistable 
system. In comparison to the previously constructed 
bistable circuit comprised of two mutually inhibiting 
switches [20], our single-switch system is simpler (four 
strands vs six) but more delicate to design. DNA hy- 
bridization reactions can no longer be made arbitrarily 
strongly favored and irreversible; the 'inhibition' and 
'release' reactions must therefore be designed carefully. 
Further, the multistep activating cascade makes use of 
both 'sense' and 'antisense' sequences, making it dif- 
ficult to eliminate spurious reactions such as 'interfer- 
ing', 'recovering', and 'recapturing'. Fortunately, the 
working solution demonstrated here ought to be easily 
generalized to other sequences, enabling the modular 
construction of in vitro transcriptional circuits contain- 
ing both activating and inhibiting switches. 

Fulfilling the promise of predictable system design, 
when configured to self-activate by using the same se- 
quence elements for both input and output domains, 
the expected bistable behavior was observed. Some- 
what more surprisingly, when a simple network model 
was trained on limited data, system sensitivities to DNA 
and enzyme concentrations were qualitatively and some- 
times quantitatively predicted. These studies also sug- 
gest that the addition of a single-stranded RNA exonu- 
cleasc, such as RNase R, would help clean up incom- 
plete degradation fragments, improve performance, and 
yield a more predictable system. 

While in principle arbitrary behaviors can be ob- 
tained using networks of inhibiting switches (for ex- 
ample, activation can be obtained by two inhibiting 
switches in series), the direct implementation of acti- 
vation using a single switch provides some advantages. 
As discussed before, it yields simpler circuits in terms 
of the number of DNA strands required. Furthermore, 
the hybridization reactions are faster and more energy 
efficient than the transcription reactions they replace. 
This observation raises a question for the future con- 
struction of larger in vitro circuits: how much of the 
desired logic should be implemented by transcriptional 
circuitry |191 120| and how much should be implemented 
as strand displacement circuitry [27l [35l [31] ? 

Materials and Methods 

DNA oligonucleotides and enzymes 

The sequence of all DNA molecules and expected 
RNA transcript sequences were chosen to minimize the 
occurrence of alternative secondary structures, checked 
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by the Vienna group's DNA and RNA folding pro- 
gram m. All DNA oligonucleotides were purchased 
from Integrated DNA Technologies (USA). In addition 
to unlabeled A and dl, A labeled with Cy5 at the 3' end 
and dl labeled with lowaBlack-RQ at the 5' end were 
purchased. The T7 RNA polymerase (enzyme mix), 
transcription buffer, and NTP were purchased as part 
of the T7 Megashortscript kit (Ambion, Austin, Texas, 
USA; #1354). RNase H (Ambion; #2293) was pur- 
chased. RNase R was a gift from Dr. Murray P. Deutscher 
at University of Miami, school of Medicine. Note that 
according to the manufacturer's patent (# 5,256,555), 
the enzyme mix of T7 Megashortscript kit contains py- 
rophosphatase to extend the lifetime of the transcrip- 
tion reaction; since pyrophosphatase is involved in reg- 
ulating the power supply for our transcriptional circuits 
and is not directly involved in the dynamics, we neglect 
this enzyme in our models and do not call it an "essen- 
tial enzyme" for the circuit dynamics. 

Transcription 

DNA templates (T-nt and T-t strands) were an- 
nealed with 10% (v/v) lOx transcription buffer from 
90°C to 20°C over 1 hour at 10 times the final concen- 
tration used. To the annealed templates, DNA activa- 
tor and inhibitor, A and dl, were added from high con- 
centration stocks (50 fjM), 7.5 mM each NTP, and 8% 
(v/v) lOx transcription buffer were added. After adding 
all ssDNA strands and RNA signals, enzymes (RNAP, 
RNase H, and RNase R) were added and mixed. Tran- 
scription reactions for spectrofluorometer experiments 
were prepared as a total volume of 60 /iL. Transcription 
reactions for gel studies were prepared as a total volume 
of 10 fiL and were stopped by denaturing dye (80% for- 
mamide, 10 mM EDTA, 0.01 g XCFF). For the purifi- 
cation of RNA signal rA for the experiments and for gel 
controls, the full length template side strand (the com- 
plement of T-nt rather than T-t) was used to prepare 
a fully duplex DNA template for rA. The transcription 
reaction was prepared as a total volume of 60 fjL with 
0.2 ^M fully duplex DNA templates. The transcription 
condition was the same as above except that no A or dl 
were added, 20% (v/v) RNAP was used, and RNases 
were omitted. After a 6 hour incubation at 37°C, the 
reaction mixture was treated with 2.5 DNase I for 
30 min to remove DNA template, and stopped by the 
addition of denaturing dye. The reaction mixture was 
run on 8% denaturing gel, RNA bands were excised and 
eluted from gel by the crush-and-soak method, ethanol 
precipitated, and rcsuspended in water. 

Data acquisition 

For spectrofluorometer experiments, excitation and 



emission for Cy5-labeled A were at 648 nm and 665 nm. 
The fluorescence was recorded every minute using a 
SPEX Fluorolog-3 (Jobin Yvon, Edison, New Jersey, 
United States) and normalized against maximum flu- 
orescence (measured in the presence of Cy5-labeled A 
prior to the addition of quencher-labeled dl) and min- 
imum fluorescence (measured after the addition of ex- 
cess quencher-labeled dl) accounting for volume increase 
due to the addition of rA and enzymes. Denaturing 
polyacrylamide gels (8% 19:1 acrylamideibis and 7 M 
urea in TBE buffer) were allowed to run for 50 min 
with 10 V/cm at 65°C in TBE buffer (100 mM Tris, 
90 mM Boric Acid, 1 mM EDTA). The 10-base DNA 
ladder (Invitrogen, Carlsbad, California, United States; 
#10821-015) was used in the control lane. The non- 
denaturing gels (10% 19:1 acrylamide:bis in TAE buffer) 
were allowed to run for 100 min with 13V/cm at 35°C 
in TAE buffer containing 12.5 mM Mg^+ (40 mM Tris- 
Acetatc, 1 mM EDTA, 12.5 mM Mg-Acetate, pH 8.3). 
The gels were stained with SYBR gold (Molecular Probes, 
Eugene. Oregon, United States; #S-11494) and the gel 
data was quantitated using the Molecular imager FX 
(Biorad, Hercules, California, United States). The to- 
tal concentrations of RNA product r A in the denaturing 
gel were measured with respect to 0.5, 1, and 1.5 yuM 
puriflcd rA in control lanes. 

Model simulation 

The kinetic simulations and parameter fittings were 
implemented in MATLAB. Differential equations were 
solved using the ode23s routine. The parameter sets 
used in Figure 6 were derived from the Metropolis' min- 
imization of sum squared differences between the exper- 
imental results in Figure 6A and the simulation results 
of kinetic model described in Appendix. We identify 
model parameters, 9, as a. vector of logs of rate con- 
stants; the kinetic model has 11 parameters {kTA, kAi, 

krAI, kTAI, kAIrA, kM,ON, kM,OFF, kcat,ON, kcat,OFF, 

kM,H^ ^iid kcat,H)- Each parameter has pre-specified 
minimum and maximum values (see Appendix). The 
cost function is defined as 

N 

where the concentration of rA is measured at 210 min in 
/iM scale both for the experiment and simulation. The 
training set (Figure 6A) contained fifty experimentally 
measured [rA] values (N = 50). Starting from random 
initial values, for each iteration, a random 11x1 vector 
of -1, 0, or -1-1 (with equal probabilities) was generated 
to decide whether to decrease, remain, or increase each 
parameter in the next step. The step size for each pa- 
rameter was 1/ 100th of the range for that parameter in 
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log scale. We accepted or rejected the next step based 
on Metropolis criteria — i.e., the acceptance probabil- 
ity was p = 1 if AE" < and p = e'^^l^ if Ai; > 
with T = 20-2 = 0.125 . Wc used three sampling 
trajectories with 620,000 iteration steps starting from 
random parameter values, discarding the first 20% of 
iterations and sampling one parameter set per 30,000 
iterations. This collection of parameter sets (M = 51) 
was used for simulation results in Figure 6; they show 
the 95% confidence interval for the parameter values — 
i.e., the four (out of 51 parameter sets) with largest E 
are omitted. See Appendix for details. 

DNA sequences 

T-nt (106 mCr), S'CATTAGTGTCGTTCGTTCACAGTAATACGACT- 

CACTATAGGGAGAAACAAAGAACGAACGACACTAATGAACTACTACTACA- 

CACTAATACTGACAAAGTCAGAAA-3\ 

T-t (79 mCr), S'TTTCTGACTTTGTCAGTATTAGTGTGTAGTAGTAGT- 

TCATTAGTGTCGTTCGTTCTTTGTTTCTCCCTATAGTGAGTCG-3'. 

A (36 mer), S'TATTACTGTGAACGAACGACACTAATGAACTACTAC-S'. 

dl (38 mer), 5"-GTGTGTAGTAGTAGTTCATTAGTGTCGTTCGTTCAC- 

AG-3-. 

rA (67 mer), s-gggagaaacaaagaacgaacgacacuaaugaacua- 

CUACUACACACUAAUACUGACAAAGUCAGAAA 3'. 

Note that the RNA activator rA sequence is iden- 
tical to II in [20], and therefore, T-t sequence is also 
identical to T12-t in [20] . 
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Appendix 



A simple kinetic model for Self- activating switch 

To mathematically explore the behavior of the self-activating switch, we developed a model based on the 
DNA and RNA hybridization reactions, and enzyme reactions. Note that for the self-activating switch, the input 
RNA signal is identical to the output RNA signal, thus both are labeled rA. We constructed a mathematical 
model without RNase R, because RNase R was not used for the results presented in Figures 5 and 6 where the 
model is used. To simplify the model, we also assumed that, although an 'interfering' reaction exists that results 
in the T-rA complex, the 'recovering' and 'recapturing' reactions involving A and dl are fast enough such that 
the concentration of the T-rA complex is negligible (see Figure 3). The hybridization and strand displacement 
reactions and enzyme reactions are summarized below. 
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Michaelis-Menten Enzyme reactions 
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We do not consider side-reactions or incomplete production and degradation products. We further simplified 
enzymatic equations from the full Michaelis-Menten equations to an approximation using the pseudo-steady- 
state assumption of enzyme-substrate complex, which is reasonably accurate when enzyme concentrations are 
low compared to substrate concentrations. 
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We express the available enzyme concentrations using the pseudo-steady-state assumption as follows: 



[RNAP] = 
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where the Michaelis constants are calculated as K 
enzymes. The effective rate constants for enzymes are as follows: 

ProdoN = ^^^ph2E. [RNAP] , ProdoFF = t'"*'"'"'' I^NAP] , Degn 



to determine the affinity of substrates to the 
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Note that this approximation procedure was also used in |20] . Taken together, a set of four ordinary differential 
equations describes the dynamics of the self-activating switch as follows. 
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A:TA[T][A]+fcTA/[T-A][dI] 

kAl[A][dl] - kTA[T][A] + kAIrA[A ■ dl] [r A] 

kAi[A][dl] - fc,A/[rA][dI] - kTAi[T ■ A][dl] + DegnirA ■ dl] 

A:,A/[rA][dI] - kAirA[A ■ dI][rA] + ProdoN[T ■ A] + ProdoFF[T] 

The remaining variables, [T-A], [A-dl], and [rA-dl] are calculated from the conservation relation as [T*°*], [A*°'], 
and [dl*°*] remain constants throughout the reaction. 

Sampling parameter space with Metropolis criteria 

Using the mathematical model described above, we chose to explore the parameter space in the vicinity of 
the best fit to data based on the Metropolis sampling method, as outlined in [6]. The experimental data shown 
in Figure 6 A is used as the training set. Other experimental data are reserved as the test set. The cost function 
is calculated as 

N 

71=1 

where the concentration of rA is measured at 210 min in /_tM scale both for the experiment and simulation. The 
training set contained fifty experimentally measured [rA] values {N = 50). To estimate errors in experimental 
measurements, most experimental conditions were measured in duplicate. For certain experimental conditions 
— e.g. in the transition region of blue curve in Figure 6A — the duplicate measurements differed as much as 
0.65 fiM. We chose cr = 0.25 /iM as a reasonable estimate of the standard deviation of repeated experimental 
measurements. 

The kinetic model has a total of 11 parameters: kxA, kAi, krAi, kxAi, kAirA, kM.ON, kM.OFF, kcat,ON, 
kcat,OFF, kM,H, and kcat,H- Wc set the lower and upper bounds for these parameters as (in logio scale) [4.5, 
4.5, 4.5, 4.5, 4.5, -8, -7, -2, -3, -7, -2] and [6.0, 6.0, 6.0, 6.0, 6.0, -6, -5, -1, -1, -5, 0], respectively In order to 
minimize the effect of these widely separated scales and avoid numerical issues, we deal with the logarithms of 
the parameters for all our calculations. Henceforth, we identify model parameters, 9, as a vector of logs of rate 
constants. We initiated the parameter set as a random 11x1 vector within these bounds. For each iteration, we 
generated a random 11x1 vector of -1, 0, or +1 (with equal probabilities) to decide whether to decrease, remain, 
or increase each parameter in the next step. The step size for each parameter was set as 1/ 100th of the range 
for that parameter in a log scale. If the simultaneous updates of 11 parameters resulted in lower cost function 
(i.e. AE < 0), the updates are accepted. On the other hand, if the simultaneous updates of parameters resulted 
in higher cost function (i.e. AE > 0), the updates are accepted with a probability of p = e~^^/^ following 
Metropolis criteria. We accepted or rejected the parameter updates based on the Metropolis criteria with T — 
0.125 (= 20-2). 

Following these sampling procedures, 620,000 iterations were performed for trajectories 1 through 3 and 
470,000 iterations were performed for trajectory 4 starting from random initial conditions. To assess whether 
there are multiple attractor basins for parameter values, we plotted parameter values for the whole trajectories 
as histograms (Figure SI) with the initial values marked by red circles. The parameter distributions for all 11 
parameters achieved similar ranges and shapes irrespective of their initial values. Therefore, we conclude that 
the Monte Carlo sampling with Metropolis criteria has converged for these trajectories. 



d[T] 
dt 

d[A] 
dt 

d[dl] 

dt 
d[rA] 

dt 
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Figure SI: Histograms for parameter-sampling trajectories. (A— D) The log-scale parameter values {6) for the whole 
sampling trajectories are plotted as histograms. The x-axes are in log scale with the respective minimum and maximum 
values corresponding to the range of each parameter as stated earlier. Following Metropolis criteria, 620,000 iterations 
were performed for trajectories 1 through 3 (A-C) and 470,000 iterations for trajectory 4 (D). The initial values for 
parameters are marked by red circles. 

Our next goal is choosing M parameter sets, in such a way that they are independent of each 

other. These parameter sets eannot be randomly chosen from the above sampling trajectories because the 
different samples may be highly correlated. To estimate decorrelation time of parameter values so as to obtain 
independent samples, we analyzed the autocorrelation function of parameter values between the 100,000th and 
200,000th iteration steps (Figure S2). For most parameters, the autocorrelation function quickly decays to zero 
— e.g., Km,on, kcat,ON, and kxAi- On the other hand, for some parameters such as Km,off and kcat,OFF, 
the parameter values did not converge and slowly drifted over time, possibly because they do not impact the 
cost function in a significant way. For these two parameters, which were poorly constrained in the histogram 
analysis, the autocorrelation function did not completely decay even after 50,000 steps. Nevertheless, the signs 
of autocorrelation function at time lags beyond 30,000 steps were roughly random — it could be positive or 
negative depending on the sampling ranges. Thus, we used 30,000 iterations as the decorrelation time. 

For sampling trajectories 1, 2, and 3 (with 620,000 iterations each), we discarded the first 20% of iterations 
and selected one parameter set per 30,000 iterations, resulting in 17 parameter sets per trajectory. This collection 
of parameter sets — 51 in total (M = 51) — was used for simulation and compared with the experimental results 
in Figure 6. The line plots in Figure 6 show the central 95% of predictions — the two parameter sets with lowest 
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Figure S2: Autocorrelation functions for parameter-sampling trajectories. (A-D) The log-scale parameter values 
[0) between the 100,000th and 200,000th iteration steps for the sampling trajectories are plotted as autocorrelation 
functions. Autocorrelation functions are for the sampling trajectory 1 (A), 2 (B), 3 (C), and 4 (D). Decorrelation 
time for each parameter can be estimated by measuring lag times beyond which the autocorrelation function is close 
to zero. 



rA output predictions and the two parameter sets with highest rA output predictions were omitted — i.e., the 
four parameter sets with largest errors. For the final 51 parameter sets, the sum squared errors ranged from 
2.03 to 2.43 fiM^ for the whole training set {N = 50 measurements). One can construct an empirical covariancc 
matrix 6 from the ensemble of parameters {Sj}m=i^ 

9 = {{9-me-{9)r), 

where the angle brackets denote ensemble average. An eigenvalue decomposition of this matrix (called principal 
component analysis (PC A)) can then be inverted and information about soft and stiff modes are obtained 
analogous to that using the Hessian because PCA Hessian P = . Mode spectrum and eigenvector projections 
are shown with eigenvalue-eigenvector correspondence indicated by the numbers 1-11 (Figure S3). 

These results show that the eigenvalue spacing is almost uniform in logspace; there is no clear cutoff be- 
tween "stiff" and "soft" eigenvalues. However, the wide spacing between eigenvalues indicate that the stiffest 
few parameter combinations can be used to explain and fit most of the variations in the experimental data. 
These "sloppy model" features have been observed in many other high-dimensional multiparameter nonlinear 
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Figure S3: Mode spectrum and eigenvector projections from PCA of inverse of covariance matrix 0. (A) Eigenvalues 
plotted in log scale. The first few eigenvalues are much larger than the rest — i.e., the few combinations of parameters 
associated with large eigenvalues can be used to explain and fit most of the variation in the experimental and simulation 
data. (B) Eigenvector projections. Eigenvector-eigenvalue correspondence is indicated by the numbers 1-11: small 
numbers correspond to large eigenvalues (stiff modes) and large numbers correspond to small eigenvalues (soft modes). 



models [6]. These characteristics are typicahy encountered in models for biochemical regulatory networks in 
biological organisms [llj . 

Let us proceed to interpret stiff modes found in the eigenvector projections. The variance in the stifFest mode 
(mode 1) was largely explained by contributions from Km,on and kcat.ON as well as Km,h and kcat.H- The 
angle on KM,ON-kcat,ON plane indicates that these two values must change in opposite directions, which would 
change the ON-state switch transcription rate significantly. At the same time, the angle on K]\j^H-kcat,H plane 
indicates that these two values also must change in opposite directions, which in this case would change the 
degradation rate significantly. Together, the signs for these four parameters in eigenvector indicates that the 
stifFest mode is captured by an increased ON-state transcription rate with a decreased degradation rate and vice 
versa — i.e., "production/degradation balance" . The second stiff mode (mode 2) was dominated by contributions 
from Km,off and kcat,OFF, indicating that the changes in opposite directions for these two parameters would 
result in significant changes in the OFF-state switch transcription rate — i.e., "leakiness". The inhibition rate 
kxAi contributed significantly in modes 3 and 4. However, since the parameters that change in the same or the 
opposite directions in relation to kxAi were not clearly identified, the eigenvalue projections involving kxAi are 
not shown in this figure. 

Now we examine some of the soft modes revealed in the eigenvector projections. The variance in the 
softest mode (mode 11) is mostly captured by parameters Km.off and kcat.OFF- The angle on K]\j^off- 
kcat,OFF plane indicates that these two values must change together to preserve the OFF-state transcription 
rate, kcat,OFF / Km,off- Similarly, the next soft mode (mode 10) is captured by parameters Km,h and kcat,Hi 
preserving the degradation rate, kcat,H / Km,h- The projections onto kxA-kAi plane revealed that the softer 
mode (mode 7) corresponds to concerted changes in the activation rate (kxA) and the annihilation rate (kAi), 
which would change the time-scale of switching but preserve the average switch states. The krAi-krAAi plane 
analogously associated the softer mode (mode 9) with the concerted changes of these two parameters. 
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